    module statistics
    implicit none
    real(8) :: sumk_mc, sumk_cmfd, sumk2_mc, sumk2_cmfd
    real(8) :: avgk_mc, avgk_cmfd, stdk_mc, stdk_cmfd, stdk_mcc
    integer nacc
contains

subroutine init_statistics
    sumk_mc = 0.
    sumk2_mc = 0.
    sumk_cmfd = 0.
    sumk2_cmfd = 0.
    avgk_mc = 0.
    avgk_cmfd = 0.
    stdk_mc = 0. 
    stdk_cmfd = 0. 
    nacc = 0;

end subroutine

subroutine add_statistics(k_mc, k_cmfd)
    real(8) :: k_mc, k_cmfd
    
    nacc = nacc + 1
    sumk_mc = sumk_mc + k_mc
    sumk2_mc = sumk2_mc + k_mc**2
    avgk_mc = sumk_mc/nacc
    sumk_cmfd = sumk_cmfd + k_cmfd
    sumk2_cmfd = sumk2_cmfd + k_cmfd**2
    avgk_cmfd = sumk_cmfd/nacc
    
    if (nacc.ne.1) then
        stdk_mc = sqrt((sumk2_mc-sumk_mc**2/nacc)/nacc/(nacc-1))
        stdk_mcc = sqrt((sumk2_mc-sumk_mc**2/nacc)/(nacc-1))
        stdk_cmfd = sqrt((sumk2_cmfd-sumk_cmfd**2/nacc)/nacc/(nacc-1))
    end if 
end subroutine    

end module